function [M,var_M] = compute_data_moments(X_monthly,Nu)
 
% Periods 1-6: pre shock
% Periods 7-14: post shock

% Replications of the bootstrap
B                                                       = Nu.B;

% Blocks used in the bootstrap

% First equation: adoption rate
X                                                       = X_monthly(7:end,:);
TimeZeroX                                               = repmat(X_monthly(1,:),[8,1]);
dX                                                      = X - TimeZeroX;
TimeShock                                               = ones(size(dX));
TimeShock(1:3,:)                                        = 0;
InitialCondition                                        = X_monthly(6,:)-X_monthly(1,:);
InitialCondition                                        = repmat(InitialCondition,[size(dX,1) 1]);
InitialCondition                                        = (InitialCondition - mean(InitialCondition(:)))/std(InitialCondition(:));

% Second equation (residuals from second equation)
Z2                                                      = [ones(size(dX(:)))];

% Third equation: cross-district adoption variance
VarXT                                                   = var(X_monthly(7:end,:)-TimeZeroX,[],2);
IndicPost                                               = ones(size(VarXT));
IndicPost(1:3)                                          = 0;
Z3                                                      = [ones(size(VarXT)) IndicPost];
Y3                                                      = VarXT;

% Fourth equation: within-district adoption variance
VarX                                                    = var(X_monthly(7:end,:)-repmat(X_monthly(1,:),[8,1]),[],1);
Z4                                                      = [ones(size(VarX(:)))];
Y4                                                      = VarX(:);

Mb                                                      = NaN(8,B);

for b=1:B
   
    % Construct residuals for third equation
    InitialConditionb                                   = InitialCondition(:,Nu.indices1(:,b));
    TimeShockb                                          = TimeShock(:,Nu.indices1(:,b));
    dXb                                                 = dX(:,Nu.indices1(:,b));
    Z1                                                  = [ones(size(dXb(:))) TimeShockb(:) InitialConditionb(:) TimeShockb(:).*InitialConditionb(:)];
    Y1                                                  = dXb(:);
    Mb1                                                 = Z1\Y1;
    resids_sq                                           = (Y1 - Z1*Mb1).^2;
    
    % Estimation of coefficients
    Yb                                                  = [Y1; resids_sq; Y3(Nu.indices3(:,b)); Y4(Nu.indices4(:,b))];
    Zb                                                  = blkdiag( Z1, Z2, Z3, Z4 );
    Mb(:,b)                                             = Zb\Yb;
    
    % Transformation of coefficients
    Mb(5,b)                                             = sqrt(abs(Mb(5,b)));
    Mb(6,b)                                             = sqrt(abs(Mb(6,b)));
    Mb(7,b)                                             = sqrt(abs(Mb(7,b)));
    Mb(8,b)                                             = sqrt(abs(Mb(8,b)));
    
end

M                                                       = mean(Mb,2);
var_M                                                   = 0;
for b=1:B
   var_M = var_M + ( (Mb(:,b)-M) )*( (Mb(:,b)-M)' );  
end
var_M   = 1/(B-1)*var_M;

end

